Parameters estimation of a noisy sinusoidal signal with time-varying 

amplitude 



Da-yan Liu, Olivier Gibaru and Wilfrid Perruquetti 



o 



< 



> 

m 
o 
cn 

in 
o 



Abstract — In this paper, we give estimators of the frequency, 
amplitude and phase of a noisy sinusoidal signal with time- 
varying amplitude by using the algebraic parametric techniques 
introduced by Fliess and Sira-Ramirez. We apply a similar 
strategy to estimate these parameters by using modulating 
functions method. The convergence of the noise error part due 
to a large class of noises is studied to show the robustness and 
the stability of these methods. We also show that the estimators 
obtained by modulating functions method are robust to "large" 
sampling period and to non zero-mean noises. 

I. INTRODUCTION 

Recent algebraic paramettic estimation techniques for 
linear systems [1], [2], [3] have been extended to various 
problems in signal processing (see, e.g., [4], [5], [6], [7], [8]). 
In [9], [10], [11], these methods are devoted to estimate the 
frequency, amplitude and phase of a noisy sinusoidal signal 
with time-invariant amplitude. Let us emphasize that these 
methods, which are algebraic and non-asymptotic, exhibit 
good robustness properties with respect to corrupting noises, 
without the need of knowing their statistical properties (see 
[12], [13] for more theoretical details). We have shown in 
[14] that the differentiation estimators proposed by algebraic 
parametric techniques can cope with a large class of noises 
for which the mean and covariance are polynomials in 
time. The robustness properties have akeady been confirmed 
by numerous computer simulations and several laboratory 
experiments. In [15], [9], modulating functions methods 
are used to estimate unknown parameters of noisy sinu- 
soidal signals. These methods have similar advantages than 
algebraic parametric techniques especially concerning the 
robustness of estimations to corrupting noises. The aim of 
this paper is to estimate the frequency, amplitude and phase 
of a noisy time-varying amplitude sinusoidal signal by using 
the previous two methods. We also show their stability by 
studying the convergence of the noise error part due to a 
large class of noises. 

In Sectionini we give some notations and useful formulae. 
In SectionHniand Section HVl we give parameters' estimators 
by using respectively algebraic parametric techniques and 
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modulating functions method. In Section |Vl the estimators 
are given in discrete case. Then, we study the influence of 
sampling period on the associated noise error part due to a 
class of noises. In Section IVll inspired by [15] a recursive 
algorithm for the frequency estimator is given, then some 
numerical simulations are given to show the efficiency and 
stability of our estimators. 

II. Notations preliminaries 



Let us denote by Dj := {T E 



[0, T] C £1}, and 



w 



li.Ki'^) = (1 -t)^t'^ for any T e [0, 1] with li,Ke] - l,+°o[. 
By using the Rodrigues formula (see [16] p. 67), we have 

dr' 

where P^^ ''^ \ min(fc,/x) > / G N, is the order Jacobi 
polynomial defined on [0,1] (see [16]): Vt £ [0,1], 



-{wm,.(t)} = (-i)'/!w^_,._,(^)^r'''"'(T), (1) 



s=0 



w,_,.,(t). (2) 



Then, we have the following lemma. 

Lemma 1: Let / be a '^"+'(£2) -continuous function (n S 
N) and 11^ ^ be a differential operator defined as follows 



1 



k.ti 



(3) 



where s is the Laplace variable, A: G N and — 1 < /i G R. 
Then, the inverse Laplace transform of 11^' ^/ where / is the 
laplace transformation of / is given by 



J^-'[lVlJ{s)]{T) 



_^n+l+pL+k 



p.+n,k 



W 



nMn^^)fiTT)dT, 



(4) 



where T eDj and c^,+„^k = Y(n+n+\] - 

In order to prove this lemma, let us recall that the a- 
order (a G M+) Riemann-Liouville integral (see [17]) of a 
real function g (K — J> R) is defined by 



1 



r(a) Jo ' 

The associated Laplace transform is given by 



(5) 



(6) 



where g denotes the Laplace transform of g. 

Proof. Let us denote W^+„,f:+„{t) = {T - tf+"t''+" for 

any r G [0,r]. Then, by applying the Laplace transform to 



the following Riemann-Liouville integral and doing some 
classical operational calculations, we obtain 



Let us apply the inverse Laplace transform to (fTZt . then by 
using Lemma [T] we obtain 

\l (w;,'|4,,+4(T)+2(a)r)2w^+4_,+4(T))x(rT) 



Then, by substituting T by Tx we have 

(T)/W(T)^T 



(7) 



+ {(OTy Wf,+4,k+^iT)x{TT)dT^0. 

According to O, we have w^^_^^i^^^{0) = m'|1+4,a.+4(1) for 
/ = 0, . . . , 3. Then by applying integration by parts, we get 



+ j^'vv^+4,i+4(T)x(4)(rTyT = 0. 

Thus, ttP' is obtained by 



By using ©, we obtain wJ;^„^+„(0) = wj;^„^^„(l) = for 
/ = O,-- - ,n — L Finally, this proof can be completed by 
applying n times integration by parts to O. □ 



III. Algebraic parametric techniques 

Let y ^ x + UJ be a noisy observation on a finite time 
interval £2 C of a real valued signal x, where CJ is an 
additive corrupting noise and 



2A.V 



(13) 



Vr G £2, x(r) = (Aq +Aif ) sin(a)f + 0) 



(8) 



with Aq eM!^, Ai eK*, 0) eM!^ and ~j7:<(j) < \n. Observe 
that X is a time-variant varying sinusoidal signal, which is a 
solution to the harmonic oscillator equation 

Vfe£2, x^'^\t)+2o?x{t) + co'^x{t)=Q. (9) -^l- /o vv^+4,i+4(T) sin(t(jrT+ 0)c/t > 0, then we get 



where = /g w^+4_(t+4(T)x(rT)t/T, 

B, = 2/Jw^+4,^+4(T)x(2)(rT) ofT, d = 

/o vv^+4,A:+4(T)x(4)(rT)t/T. Since x('*)(rT)+2coV2)(rT) + 
(B'*x(rT) — for any T G [0, 1], we get 

1(b2_4A,4) = 

/■I \ ' 

W^+4,i:+4(T)x(2)(rT) dX + (O^W ^+4^k+A{T:)x{T x) dx\ 

Observe that xP) [Tx] + op-x{Tx) ^ 2coAi cos{coTt + <j)) for 
any xe [0,1]. If «Ai /J w^+4,i+4(T) cos(a)rf 



Then, we can estimate the parameters co, Aq and (j) by 
applying algebraic parametric techniques to (|9}. 

Proposition 1: Let k G N, -1 < /i e K and T e Dj 
such that Ai vv'^+4jt+4(T) sin(a)rT + 0)£/t < 0, then the 
parameter o is estimated from the noisy observation y by 



-B,+ ^Bl~AA,C, 



0) 



(14) 



0) 



By+JB^-AAyCy 



2Av 



(10) 



Finally, this proof can be completed by applying integration 
by parts and substituting x by y in the last equation. □ 
By observing that xq = x(0) = Aosin^, xq = 
i(0) — AqOJcos^ + Aisin^ and Xq^' = x('')(0) = 
—O^XQ — 2a)^Aisin0, then we can obtain Aqcos^ = 



where A,. = jQW^+4^t+A{T:)y{Tx)dx, By = 



2T^J^w^+4M4{^)yiT^)dx, Cy = !^w'^l,,^,{x)y{Tx)dx, 



1 ,.,(4) 



is given by ([TJ with / = 1,2,4. 



Proof. By applying the Laplace transform to (|9]l, we get 

/x(.s) + 2co^s'^x{s) + co'^xis) 



=s^xo + s^XQ + {2co^xo+xq)s + {2co^XQ +xl^^) 



(11) 




= arctan 



Let us apply A; + 4 (A; e N) times derivations to both sides of 
(fTTT i with respect to s. By multiplying the resulting equation 
by s^^^^ with -1 < jU e M, we get 



n2,^i(^) + 2a)2n2+2^^+2x(^) + a)*n^V4.M+4^^(^) = «■ d^) 



Thus, we need to estimate xq, xq and Xq so as to obtain the 
estimations of Aq and (p. 

Proposition 2: Let — 1 < /X G K and T G Dj, then the pa- 
rameters Ao and are estimated from the noisy observation 



y and the estimated value of co given in ( fTOl i: 



Ao=\xl 



2\ I 



where 




XQ^ [ 




Jo 




-(3) 1 






/' 

Jo 
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-5) 



4«6 

= arctan — 

T Jo 



-4(^n'I0)a±i±a_i»„««3-,(.) 

- (a)r)'^C^+43W^+4,3(T), 
-2 

■Pf (t) = c^,3W^3(t) + llc^+1, 2^^+1,2(7) 



(16) 



r(Ai + 6) 

+ 28c^+2aW^+2,l(T) + 12c^+3_OW^+3,o(t) 
+ 2{a)Tf (c^+2,3W^+2,3(t) +5c^+3,2W^+3,2(t)) 
+ 4(d)r)^ (C|j+4,1H';i+4,i(t) ~C^+5,OW^+5,o(t)) 
+ (aT)'^ (c^+43W|i+43(T) - C^+5,2W|i+5,2(T)) , 



-6 



^ ^3\3!c„+,3- 



i=0 



1=0 



Proof. In order to estimate xq, we apply the following 

1 ^ 
operator lli ~ — -p- • — to (lllb with -1 < u £ M, which 

annihilates each terms containing x'q for / = 1,2,3. Then, by 
using the Leibniz formula, we get 



3\ 4! 1 



■3\ 2! 1 



f(3-0 



0) 



Let us express the last equation in the time domain. By 
denoting T as the length of the estimation time window we 
have 



r(M+5) io ,eoV'7 (4-0 



+(0^^ C^i+4,3W^+4,3iT:)dT. 



Hence, by substituting T by Tt, x by y and taking the esti- 
mation of (0 given in Proposition [T] we obtain an estimate for 

■, , , , ^ 1 d I d^ 
xq. Similarly, we apply the operator 112 — — rj • — rn 

5^'+^ ds s ds^ 

1 d I 

(resp. 113 = — —7 • — T • -) to (fTTI) to compute an estimate 

s^^+* ds^ s 

(3) 

for io (resp. x\ ). Finally, we get estimations for Aq and ^ 
from relations (fTsT i by using the estimations of xq, xq, Xq 
and (0. □ 

IV. Modulating functions method 

Proposition 3: Let / be a function belonging to '^'*([0, 1]) 
which satisfies the following conditions /'''(O) = /'''(I) for 
i = 0, . . . , 3. Assume that Ai /J /(t) sin(orT + 0)iiT < with 
T G ZJj-, then the parameter (O is estimated from the noisy 
observation y by 



CO = 



-By+^Bj-AAyCy 



(17) 



where Ay = f{x)y{Tx)dx, B, = 2T\f^l f{T)y{TT)dt, 
Cy = S^f('\x)y{Tx)dT. 

Proof. Recall that x^^^ {Tx)+2(0^x^^^ (Tt) + co'^x{Tt) = for 
any T G [0, 1]. As / is continuous on [0, 1], then we have 

/' /(t)x(^' {TT)dT + 20)2 t f(T)x^^^ {Tx)dx 
Jo Jo 

+ 0)'* / f{T)x{Tx)dT = 0. 
Jo 

Then, this proof can be completed similarly to the one of 
Proposition [T] □ 
Proposition 4: Let / for / = 1 , . . . , 4 be four continuous 
functions defined on [0, 1]. Assume that there exists T G Dj 
such that the determinant of the matrix = (^f})i</j<4 
is different to zero, where for / = 1 , ... ,4 

Mf'x = /'/i(T)sin(«rT)c/T,Mf3 = fi{T)Ttsm{(0TT)dz, 
Jo ' Jo ' 

Mi2 = f MT:)cOs{C0TT)dT,Mf^^ [ fi{T)TTCOs{(0TT)dT. 

Jo ' Jo 

Then, for any g] — f , f [ the estimations of Aq, Ai and (j) 
are given by 



A,= 

^ = arctan 



A,cos^) +(A,sin0)' 
/ Aflsin^ 



(18) 



\Aocos0 ^ 

where the estimates of A, cos0 and A, sin0 for / = 0, 1 are 
obtained by solving the following linear system 



( Aqcos^ ^ 








Ao sin 






A\ COS0 








Y Ai sin0 1 






) 



(19) 



where ij, = Jq fi{x)y{Tx)dx for i — 1,...,4, and d) is the 
estimate of (O given by Proposition [3] 



/ AqCOS^ ^ 






A() sin (p 




hi 


Ai cos(p 






\ Ai sin0 J 







Proof. Let us take an expansion of x 

x{Tt) —Aqcos^ sin(a)rT) +Aosin0 cos(a)rT) 

+A\ cos rTsin(a)rT) +Ai sm(j)TTcos{a)TT), 

where T € [0, 1], T G Dt- By multiplying both sides of the 
last equation by the continuous functions for ; = 1,...,4 
and by integrating the resulting equations between and 1, 
we obtain 

Ij. = Ao cos 0Mfi +Ao sin ^Mfj + A i cos ^Mg + A i sin 0Mf4 . 
Then, it yields the following hnear system 



Mr, 



Since det(M(u) ^ 0, we obtain A; cos (j) and A,- sin for / — 
0, 1. Finally, the proof can be completed by substituting x by 
y in the so obtained formulae of A, cos0 and A, sin0. □ 

From now on, we choose functions j^^,, with n G 

N, jUjif g] — for the previous modulating functions. 

Consequently, the estimate for co given in Proposition |3] 
generahzes the estimate given in Proposition [T] 

V. Analysis of the errors due to the noise and 

THE sampling PERIOD 

A. Two different sources of errors 

Let us assume now that y{ti) = x(f,) + £iJ(fi) (f; G i2) is a 
noisy measurement of x in discrete case with an equidistant 
sampling period Tj. Since y is a discrete measurement, 
we apply the trapezoidal numerical integration method to 
approximate the integrals used in the previous estimators. Let 
T,- = ^ and fl; > for i — 0, . . . ,m with m = ^ G N* (except 
for flo > and «,„ > 0) be respectively the abscissas and the 
weights for a given numerical integration method. Weight 
(resp. fl„,) is set to zero in order to avoid the infinite value 
at T = when - 1 < JC < (resp. T = 1 when - 1 < < 0). 
Let us denote by q the functions obtained in the integrals of 
our estimators. Then, we denote by Iq := Jq q{T)y{TT)dT. 

a- 

Hence, is approximated by 7^'™ := ^ — q{Ti)y{TTi). By 

1=0 

Iq +eg , where 



writing y(f,) =x(f,) + Cl(f;), we get 
, m 



^m.m _ ^ — q[ti)(U{Tti). Thus the integral 7;^ is corrupted 

1=0 '"" 

by two sources of errors: 

• the numerical error which comes from the numerical 
integration method, 

• the noise error contributions e^''". 

In the next subsection, we study the choice for the sampling 
period so as to reduce the noise error contributions. 

B. Analysis of the noise error for different stochastic pro- 
cesses 

We assume in this section that the additive corruption noise 
{CJ(f,),r,- G £1} is a continuous stochastic process satisfying 
the following conditions 



(Ci) : for any i,r > 0, s^ t, (U{s) and CJ(f) are indepen- 
dent; 

(C2) : the mean value function of {CJ(t),t > 0} belongs 
to 

(C3) : the variance function of {CT(t),t > 0} is bounded 
on Q.. 

Note that white Gaussian noise and Poisson noise satisfy 
these conditions. When the value of T is set, then — ?■ 
is equivalent to m +00. We are going to show the 
convergence of the noise error contributions when Ts 0. 

Lemma 2: Let n7(f,) be a sequence of {?(J(t),t > 0} with 
an equidistant sampling period 7^, where {(D{t),T > 0} be 
a continuous stochastic process satisfying conditions (Ci ) — 
(C3). Assume that q G ^^([0, 1]), then we have 



lim / q{T)E[a}{TT)]dT, 

Hm Var le^''"] = 0. 

Proof. Since n7(f,) is a sequence of independent random 
variables (Ci), then by using the properties of mean value 
and variance functions we have 



- (=0 



1 m 

Var [el'-] =—Y,aJq\Ti)Var[w{Tx,)] . 



(21) 



1=0 



According to (C3), the variance function of U5 is bounded. 
Then we have 

< -,t^hH^i) \Var[U5{TT,)]\ < t'^'i^^^ 
m^ " m " m 

(22) 

where a{m) = max a, and U = sup jVar [Cl(rT)] | < +00. 

0<Km 0<T<1 

Moreover, since q G 

^2([o,l]) and the mean value function 
of CJ is integrable (C2), then we have 

lini E[el'"] ^ (\{x)E[U5{Tx)]dT, 
lim Y-q^iTi)^ q^{T)dx<+oo. 



(23) 



As all fl, are bounded, we have JJ^^ Y —q^(Xi) —0. This 

proof is completed. □ 
Theorem 1: With the same conditions given in Lemma |2] 
we have the following convergence 

e^-'"^'-^^^^ [\(T)E\(D(TT)]dT, whenr,^0. (24) 
^ Jo 

Moreover, if noise 07 satisfies the following condition 

n-l 

(C4) : E[w{x)] = ^ V,-t' with n G N and v,- G K, 

1=0 

and q = w|j"|„ ^r+n with jj.,K e] — 5, +°°[, then we have 

lim E [ef"] = 0, (25) 

and 



0, when 0. 



(26) 



Proof. Recall that E {Y„, - cf = Var [F„,] + {E [Y,„] - cf 
for any sequence of random variables F,,, with c G R, then 
by using Lemma |2l e^'"' converges in mean square to 

n-l 

jQq{x)E[w{Tx)]dT when T, ^0. If E[w{x)] = ^V,'T' 

i=0 

and II, K g] — then by using the Rodrigues for- 

mula given by ([TJ we obtain Wu^„ xr+n 

G=Sf2([0,l]) and 

So'^^ji+n.K+n{'^)E['^[Tx)\dT — 0. Hence, this proof is com- 
pleted. □ 

VI. Numerical implementations 

In our identification procedure, we use a moving inte- 
gration window. Hence, the estimate of o at f, is given by 
Proposition |3] as follows 



2Ay,. 2A 



0,1, 



(27) 



where A,,,_. = ^B?,^ -4A„_.C,,. A,,_. = T^/}''^ , B,,^ = 



/ 



/^;:,7 with y,, 



y(r • +f,). Note that if 



Ay J, = 0, then there is a singular value in p7|) . If we denote by 

' Ay 

0, = ^— where Dy, = —By, or Dy, = Ay, , then we can apply 

}tj ' ' ' ' 

the following criterion (see [15]) to improve the estimation 
of O) 



(28) 



where / = 0, 1, . . . , and v g]0, 1]. The parameter v represents 
a forgetting factor to exponentially discard the "old" data 
in the recursive schema. The value of 0;, which minimizes 



the criterion (|28]) . is obtained by seeking the value which 
cancels Thus, we get 



0, = - 



£v'+l-^Dy,Ay,_. 

7=0 



(29) 



;=o 



'■+1-./ 



Similarly to [15], we can get the following recursive algo- 
rithm for (|29l) 



0, 



a,0,+Dy, Ay, 



0,1,. 



(30) 



where a,- = ^ v'^^^-' (Ay,, j Moreover, a,+ i can be recur- 

7=0 

sively calculated as follows a,+ i = v ( a, + ^Ay, 

Example 1: According to Section [V] we can reduce the 
noise error part in our estimations by decreasing the sampling 
period. Hence, let (y(f,) = x(f,) + cCJ(f,)),>Q be a generated 
noise data set with a small sampling period Ts = 5n x 10^"* 
in the interval [0;3;r] (see Fig. [T]i where 



sin(10f, + f), 
4t.) = { |sin(10f, + f), 

2sin(10f, + f), if2n<ti<3n, 



if < f, < n, 
if TT < f; < 27r, (31) 



and noise c{n{xi) is simulated from a zero-mean white Gaus- 
sian iid sequence with c = 0.1. Hence, the signal-to-noise 

is equal to SNR = 20.8dB. 



ratio 5A^/?=101og,o(jS|^) 
In order to estimate the frequency, by applying the previous 
recursive algorithm we use Proposition [T] with K = jj. = Q, 
m =450 and v = 1. The relating estimation error is shown in 
Fig. [21 By using the estimated frequency value, we estimate 
the amplitude and phase of the signal by applying Proposition 
|2] with H =0, m = 500 and Proposition |4] with m = 500, 
/i = W3,2, fi = W2,3, /a = W3,4 and /4 = W4,3. The relating 
estimation errors are shown in Fig. |3] and Fig. 2] We can 
observe that with small value of the relating estimation 
errors are also small. 

Example 2: In this example, we increase the value of 
to Ts — 2n X 10^^ and reduce the noise level to c = 0.01. 
Moreover, we add a bias term perturbation ^ = 0.25 in 
( I3TI 1 when f, G]27r,37i;]. The estimations of (O are obtained 
by Proposition [T] with K = ji = 0, m = 12 and v = 1. The 
estimations of the amplitude and phase are given by applying 
Proposition |2] with jU = 0, m = 12 and Proposition |4] with 

m= 15, /i = Wj']' /2 = ^^2*3' /3 = ^3 4 /4 = ^4^3' ^he 

relating estimation errors are shown in Fig. |5] and Fig. 
We can observe that the estimators obtained by modulating 
functions method are more robust to the sampling period 
and to the non zero-mean noise than the ones obtained by 
algebraic parametric techniques. 




Fig. 1. The noisy observation y and the signal x 
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Fig. 2. Relating estimation error of a 
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Fig. 3. Relating estimation en'ors of Aq 
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Fig. 4. Relating estimation errors of i 
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Fig. 5. Relating estimation en'ors of Aq 
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Fig. 6. Relating estimation errors of 



VII. CONCLUSIONS AND FUTURE WORKS 

In this paper, two methods are given to estimate the fre- 
quency, amplitude and phase of a noisy sinusoidal signal with 
time-varying amplitude, where the estimates are obtained 
by using integrals. There are two types of errors for these 
estimates: the numerical error and the noise error part. Then, 
the convergence in mean square of the noise error part is 
studied. A recursive algorithm for frequency estimator is 
given. In numerical examples, we show some comparisons 
between the two proposed methods. Moreover, these methods 
can also be used to estimate the frequencies, the amplitudes 
and the phases of two sinusoidal signals from their noisy 
sum (see [11]). The analysis for colored noises will be done 
in a future work. 
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